/**
 * @file Intp_Solver.h
 * @brief An Interpolation Solver for Multigrid Solver.
 * It can support other multigrid methods.
 * @author Linhuo (1066117119@qq.com)
 * @version 1.0
 * @date 2021-06-09
 * 
 * @copyright Copyright (c) 2021  Linhuo
 */
#ifndef _PRO2_INTP_SOLVER_H_
#define _PRO2_INTP_SOLVER_H_
#include <iostream>
#include <eigen3/Eigen/Core>
#include <cmath>
typedef Eigen::Matrix<double, Eigen::Dynamic, 1> VectorXd;

class Intp_Solver
{
public:
    virtual VectorXd Interpolate(VectorXd _r) = 0;
};

class Intp_1D_Linear : public Intp_Solver
{
public:
    VectorXd Interpolate(VectorXd _r)
    {
        // std::cout << "Linear1D" << std::endl;
        VectorXd _u_temp;
        int N = (_r.rows() - 1) * 2;
        _u_temp.setZero(N + 1);
        for (int i = 1; i <= N - 1; i = i + 2)
            _u_temp(i) = (_r((i - 1) / 2) + _r((i + 1) / 2)) / 2.0;
        for (int j = 0; j <= N; j = j + 2)
            _u_temp(j) = _r(j / 2);
        _u_temp(0) = _r(0);
        _u_temp(N) = _r(N / 2);
        return _u_temp;
    };
};

class Intp_1D_Quadratic : public Intp_Solver
{
public:
    VectorXd Interpolate(VectorXd _r)
    {
        // std::cout << "1D_Quadratic" << std::endl;
        VectorXd _u_temp;
        int N = (_r.rows() - 1) * 2;
        _u_temp.setZero(N + 1);
        _u_temp(1) = _r(0) / 2.0 + _r(1) / 2.0;
        for (int i = 3; i <= N - 3; i = i + 2)
        {
            _u_temp(i) = 9.0 * _r((i - 1) / 2) / 16.0 + 9.0 * _r((i + 1) / 2) / 16.0 - 1.0 * _r((i + 3) / 2) / 16.0 - 1.0 * _r((i - 3) / 2) / 16.0;
        };
        _u_temp(N - 1) = _r(N / 2 - 1) / 2.0 + _r(N / 2) / 2.0;
        for (int j = 0; j <= N; j = j + 2)
            _u_temp(j) = _r(j / 2);
        _u_temp(0) = _r(0);
        _u_temp(N) = _r(N / 2);
        return _u_temp;
    };
};

class Intp_2D_Linear : public Intp_Solver
{
public:
    VectorXd Interpolate(VectorXd _r)
    {
        // std::cout << "Linear2D" << std::endl;
        VectorXd _u_temp;
        int N = (sqrt(_r.rows()) - 1) * 2;
        _u_temp.setZero((N + 1) * (N + 1));
        for (int j = 1; j <= N - 1; j++)
        {
            for (int i = 1; i <= N - 1; i++)
            {
                if (i % 2 == 1 && j % 2 == 1)
                {
                    _u_temp(j * (N + 1) + i) = 1.0 / 4 * (_r((j - 1) / 2 * (N / 2 + 1) + (i - 1) / 2) + _r((j - 1) / 2 * (N / 2 + 1) + (i + 1) / 2) + _r((j + 1) / 2 * (N / 2 + 1) + (i - 1) / 2) + _r((j + 1) / 2 * (N / 2 + 1) + (i + 1) / 2));
                }
                else if (i % 2 == 0 && j % 2 == 0)
                {
                    _u_temp(j * (N + 1) + i) = _r(j / 2 * (N / 2 + 1) + i / 2);
                }
                else if (i % 2 == 0 && j % 2 == 1)
                {
                    _u_temp(j * (N + 1) + i) = 1.0 / 2 * (_r((j + 1) / 2 * (N / 2 + 1) + i / 2) + _r((j - 1) / 2 * (N / 2 + 1) + i / 2));
                }
                else if (i % 2 == 1 && j % 2 == 0)
                {
                    _u_temp(j * (N + 1) + i) = 1.0 / 2 * (_r((j) / 2 * (N / 2 + 1) + (i - 1) / 2) + _r((j) / 2 * (N / 2 + 1) + (i + 1) / 2));
                };
            };
        };

        return _u_temp;
    };
};

class Intp_2D_Quadratic : public Intp_Solver
{
public:
    VectorXd Interpolate(VectorXd _r)
    {
        // std::cout << "2D_Quadratic" << std::endl;
        VectorXd _u_temp;
        int N = (sqrt(_r.rows()) - 1) * 2;
        _u_temp.setZero((N + 1) * (N + 1));
        for (int j = 2; j <= N - 2; j++)
        {
            for (int i = 2; i <= N - 2; i++)
            {
                if (i % 2 == 0 && j % 2 == 0)
                {
                    _u_temp(j * (N + 1) + i) = _r(j / 2 * (N / 2 + 1) + i / 2);
                }
                else if (i % 2 == 0 && j % 2 == 1)
                {
                    _u_temp(j * (N + 1) + i) =
                        9.0 / 16 * (_r((j + 1) / 2 * (N / 2 + 1) + i / 2) + _r((j - 1) / 2 * (N / 2 + 1) + i / 2)) -
                        1.0 / 16 * (_r((j + 3) / 2 * (N / 2 + 1) + i / 2) + _r((j - 3) / 2 * (N / 2 + 1) + i / 2));
                }
                else if (j % 2 == 0 && i % 2 == 1)
                {
                    _u_temp(j * (N + 1) + i) =
                        9.0 / 16 * (_r((j) / 2 * (N / 2 + 1) + (i - 1) / 2) + _r((j) / 2 * (N / 2 + 1) + (i + 1) / 2)) -
                        1.0 / 16 * (_r((j) / 2 * (N / 2 + 1) + (i - 3) / 2) + _r((j) / 2 * (N / 2 + 1) + (i + 3) / 2));
                }
                else if (j % 2 == 1 && i % 2 == 1)
                {
                    _u_temp(j * (N + 1) + i) =
                        10.0 / 32 * (_r((j - 1) / 2 * (N / 2 + 1) + (i - 1) / 2) + _r((j - 1) / 2 * (N / 2 + 1) + (i + 1) / 2)) - 1.0 / 32 * (_r((j - 1) / 2 * (N / 2 + 1) + (i - 3) / 2) + _r((j - 1) / 2 * (N / 2 + 1) + (i + 3) / 2)) +
                        10.0 / 32 * (_r((j + 1) / 2 * (N / 2 + 1) + (i - 1) / 2) + _r((j + 1) / 2 * (N / 2 + 1) + (i + 1) / 2)) - 1.0 / 32 * (_r((j + 1) / 2 * (N / 2 + 1) + (i - 3) / 2) + _r((j + 1) / 2 * (N / 2 + 1) + (i + 3) / 2)) -
                        1.0 / 32 * (_r((j - 3) / 2 * (N / 2 + 1) + (i - 1) / 2) + _r((j - 3) / 2 * (N / 2 + 1) + (i + 1) / 2)) - 1.0 / 32 * (_r((j + 3) / 2 * (N / 2 + 1) + (i - 1) / 2) + _r((j + 3) / 2 * (N / 2 + 1) + (i + 1) / 2));
                    ;
                };
            };
        };

        for (int j = 2; j <= (N - 2); j += 2)
        {
            for (int i = 1; i <= (N - 1); i += (N - 2))
            {
                _u_temp(j * (N + 1) + i) = 1.0 / 2 * (_r((j) / 2 * (N / 2 + 1) + (i - 1) / 2) + _r((j) / 2 * (N / 2 + 1) + (i + 1) / 2));
            };
        };

        for (int j = 1; j <= (N - 1); j += (N - 2))
        {
            for (int i = 2; i <= (N - 2); i += 2)
            {
                _u_temp(j * (N + 1) + i) = 1.0 / 2 * (_r((j + 1) / 2 * (N / 2 + 1) + i / 2) + _r((j - 1) / 2 * (N / 2 + 1) + i / 2));
            };
        };
        for (int j = 1; j <= (N - 1); j += 2)
        {
            for (int i = 1; i <= (N - 1); i += (N - 2))
            {
                _u_temp(j * (N + 1) + i) = 1.0 / 4 * (_r((j - 1) / 2 * (N / 2 + 1) + (i - 1) / 2) + _r((j - 1) / 2 * (N / 2 + 1) + (i + 1) / 2) + _r((j + 1) / 2 * (N / 2 + 1) + (i - 1) / 2) + _r((j + 1) / 2 * (N / 2 + 1) + (i + 1) / 2));
            };
        };
        for (int j = 1; j <= (N - 1); j += (N - 2))
        {
            for (int i = 1; i <= (N - 1); i += 2)
            {
                _u_temp(j * (N + 1) + i) = 1.0 / 4 * (_r((j - 1) / 2 * (N / 2 + 1) + (i - 1) / 2) + _r((j - 1) / 2 * (N / 2 + 1) + (i + 1) / 2) + _r((j + 1) / 2 * (N / 2 + 1) + (i - 1) / 2) + _r((j + 1) / 2 * (N / 2 + 1) + (i + 1) / 2));
            };
        };
        return _u_temp;
    };
};

#else
//DO NOTHING
#endif